home *** CD-ROM | disk | FTP | other *** search
/ Amiga Tools 5 / Amiga Tools 5.iso / grafik / converter / limbo4.0 / src / 3d / findpool.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-07-16  |  7.4 KB  |  256 lines
  1. #include "includes.h"
  2.  
  3.  BitMap3D *FindPool_SampleImg3D;
  4.  
  5.  
  6.  int GridQuant(float Dist, int Size)
  7.   {
  8.    float s=(float)(Size)/2.0;
  9.    int i=(int)(Dist+s-0.5);
  10.    if (i>=Size) i=Size-1;
  11.    if (i<0)     i=0;
  12.    
  13.    return i;
  14.   }
  15.  
  16. /**************************************|****************************************
  17. Routine   :FindNodes3D
  18. Input  par:BitMap3D *Image 
  19.            PoolStructure *Pool,int B,int delta,
  20.                           int FeatureDim,float fTShadeEdge, int domain
  21. Output par:
  22. Function  :
  23. ***************************************|***************************************/
  24.  
  25.  PoolNode ****FindNodes3D(BitMap3D *Image,PoolStructure *Pool,int B,int delta,
  26.                           int FeatureDim,float fTShadeEdge, int domain)
  27.   {
  28.    PoolNode ****Array;
  29.    PoolNode *Node;
  30.    float *Features;
  31.    int x,y,z,XMax,YMax,ZMax,Count=0;
  32.    
  33.    if (domain) XMax=((Image->XSize<<1)-B)/delta +1;
  34.    else        XMax=(Image->XSize)/B;
  35.    if (domain) YMax=((Image->YSize<<1)-B)/delta +1;
  36.    else        YMax=(Image->YSize)/B;
  37.    if (domain) ZMax=((Image->ZSize<<1)-B)/delta +1;
  38.    else        ZMax=(Image->ZSize)/B;
  39.    
  40.    Array=GimmeAPoolNodeArray3D(XMax,YMax,ZMax);
  41.    
  42.    if (domain)
  43.     {
  44.      Pool->DomainX=XMax;
  45.      Pool->DomainY=YMax;
  46.      Pool->DomainZ=ZMax;
  47.     }
  48.    else
  49.     {
  50.      Pool->RangeX=XMax;
  51.      Pool->RangeY=YMax;
  52.      Pool->RangeZ=ZMax;
  53.     }
  54.    
  55.    for(x=0;x<XMax;x++)
  56.    for(y=0;y<YMax;y++)
  57.    for(z=0;z<ZMax;z++)
  58.     {
  59.      Array[x][y][z]=Node=GimmeAPoolNode();
  60.      
  61.      if (domain) Features=Classify3D(Image,(x*delta)>>1,(y*delta)>>1,
  62.                                      (z*delta)>>1,B>>1,FeatureDim,fTShadeEdge);
  63.      else        Features=Classify3D(Image,x*delta,y*delta,
  64.                                      z*delta,B,FeatureDim,fTShadeEdge);
  65.  
  66.      Node->Features=Features;
  67.      Node->Distance=NULL;
  68.      
  69.      if (Features[VAR]>fTShadeEdge)  
  70.       {
  71.        Node->x=x*delta;
  72.        Node->y=y*delta;
  73.        Node->z=z*delta;
  74.        Count++;
  75.       }
  76.      else /* decode grey range block to original domain image */
  77.      if(!domain)
  78.       {
  79.        register unsigned int X,Y,Z;
  80.        register int Par=Features[MEAN];
  81.        
  82.        for(X=x*delta;X<x*delta+B;X++)      
  83.        for(Y=y*delta;Y<y*delta+B;Y++)        
  84.        for(Z=z*delta;Z<z*delta+B;Z++) Image->Map[X][Y][Z]=Par;
  85.       }
  86.     }
  87.     
  88.    vprintf(stderr,"block size %3d, edges: %6d total: %6d  = %3.1f%c ",B,Count,
  89.                    XMax*YMax*ZMax,Count*100.0/(XMax*YMax*ZMax),'%');
  90.    
  91.    return Array;
  92.   }
  93.  
  94.  FeatureSpace *FindSpace3D(PoolStructure *Pool,int FeatureDim,int Size,float fTShadeEdge)  
  95.   {
  96.    FeatureSpace *S=GimmeAFeatureSpace(FeatureDim,Size);
  97.    Grid         *Tmp;
  98.    GridTerm     *T;
  99.    PoolNode ****DPool=Pool->DomainNodes;
  100.    PoolNode ****RPool=Pool->RangeNodes;
  101.    PoolNode    *Temp;
  102.    ListNode *TempNode;
  103.    
  104.    float *Count=GimmeAFloatArray(FeatureDim);
  105.    float *Mean =GimmeAFloatArray(FeatureDim);
  106.    float *Var  =GimmeAFloatArray(FeatureDim);
  107.    float *InvDev=GimmeAFloatArray(FeatureDim);
  108.    
  109.    float temp;
  110.    double t;
  111.    register int x,y,z,i,itemp;
  112.    
  113.    for(i=0;i<FeatureDim;i++)   /* find mean for all (domain) features */
  114.     {
  115.      double mean=0.0,var=0.0;  /* find variance cov_ii for all features */
  116.      Count[i]=0.0;
  117.      for(x=0;x<Pool->DomainX;x++)
  118.      for(y=0;y<Pool->DomainY;y++)
  119.      for(z=0;z<Pool->DomainZ;z++)
  120.       {
  121.        Temp=DPool[x][y][z];
  122.        if (Temp!=NULL && Temp->Features[VAR]>fTShadeEdge)
  123.         {
  124.          t=(double)Temp->Features[i+PREDEFFEATURES];
  125.          mean += t;
  126.          var += t*t;
  127.          Count[i]++;
  128.         }
  129.       }
  130.      Mean[i]=(float)(mean/(double)Count[i]);
  131.      t=(var-Count[i]*Mean[i]*Mean[i])/(Count[i]-1.0);
  132.      Var[i]=(float)t;
  133.      InvDev[i]=1.0/(float)sqrt(t);
  134.     }
  135.    
  136.    for(x=0;x<Pool->DomainX;x++) /* insert domain blocks in feature grid */
  137.    for(y=0;y<Pool->DomainY;y++) /* using Mahalanobis distance    */
  138.    for(z=0;z<Pool->DomainZ;z++) 
  139.     {
  140.      Temp=DPool[x][y][z];
  141.      if (Temp!=NULL && Temp->Features[VAR]>fTShadeEdge)
  142.       {
  143.        Temp->Distance=GimmeALongArray(FeatureDim);
  144.        Temp->Index   =GimmeAIntArray(FeatureDim);
  145.        Tmp=S->GridArray;
  146.        for(i=0;i<FeatureDim;i++)
  147.         {
  148.          temp=(Temp->Features[i+PREDEFFEATURES]-Mean[i])*InvDev[i];
  149.          Temp->Distance[i]=(long)(temp*F2LONG);
  150.          itemp=GridQuant(temp,Size);
  151.          Temp->Index[i] = itemp;
  152.          Tmp=Tmp->Next[itemp];
  153.         }
  154.        
  155.        T=(GridTerm *)Tmp;
  156.        if (T->Class==NULL)
  157.         {
  158.          T->Class=GimmeAListNode();
  159.          TempNode=T->Class;
  160.         }
  161.        else
  162.         {
  163.          TempNode=T->Class;
  164.          while (TempNode->Next!=NULL) TempNode=TempNode->Next;
  165.          
  166.          TempNode->Next=GimmeAListNode();
  167.          TempNode=TempNode->Next;
  168.         }
  169.        
  170.        TempNode->Domain=Temp;
  171.        T->Count++;
  172.       }
  173.     }
  174.    
  175.    for(x=0;x<Pool->RangeX;x++) /* insert range blocks in feature grid */
  176.    for(y=0;y<Pool->RangeY;y++) /* using Mahalanobis distance    */
  177.    for(z=0;z<Pool->RangeZ;z++) 
  178.     {
  179.      Temp=RPool[x][y][z];
  180.      if (Temp->Features[VAR]>fTShadeEdge)
  181.       {
  182.        Temp->Distance=GimmeALongArray(FeatureDim);
  183.        Temp->Index   =GimmeAIntArray(FeatureDim);
  184.        for(i=0;i<FeatureDim;i++)
  185.         {
  186.          temp=(Temp->Features[i+PREDEFFEATURES]-Mean[i])*InvDev[i];
  187.          Temp->Distance[i]=(long)(temp*F2LONG);
  188.          itemp=GridQuant(temp,Size);
  189.          Temp->Index[i]=itemp;
  190.         }
  191.       }
  192.     }
  193.    
  194.    FreeMeAFloatArray(Count);
  195.    FreeMeAFloatArray(Mean);
  196.    FreeMeAFloatArray(Var);
  197.    FreeMeAFloatArray(InvDev);
  198.    
  199.    return S;
  200.   }
  201.  
  202.  
  203. /**************************************|****************************************
  204. Routine   : *FindPool
  205. Input  par: BitMap *Image, int D, int B, int delta, int keepdelta
  206. Output par: none
  207. Function  : Finds the pool structure, calls FindNodes recursively finding domain
  208.             and range pool for every resolution 
  209. ***************************************|***************************************/
  210.  
  211.  PoolStructure *FindPool3D(BitMap3D *Image,int FeatureDimensions, int GridRes,
  212.                            int TShadeEdge, int NSquare,int B,int delta)
  213.   {
  214.    if (NSquare>=0)
  215.     {
  216.      PoolStructure *Pool=GimmeAPoolStructure();
  217.      register unsigned int x,y,z,i,j,k;
  218.      register int Sum;
  219.      static int first=TRUE;
  220.      
  221.      vprintf(stderr,"\n   Range  ");
  222.      Pool->RangeNodes=FindNodes3D(Image,Pool,B,B,FeatureDimensions,(float)TShadeEdge,FALSE);
  223.      
  224.      Pool->Next=FindPool3D(Image,FeatureDimensions,GridRes,
  225.                            TShadeEdge,NSquare-1,B>>1,delta);
  226.  
  227.      if(first) /* build a sampled version of the original image */
  228.       {
  229.        first=FALSE;
  230.        FindPool_SampleImg3D=GimmeABitMap3D(Image->XSize>>1,Image->YSize>>1,
  231.                                            Image->ZSize>>1,Image->ImgType);
  232.        for(x=0;x<Image->XSize;x+=2)      
  233.        for(y=0;y<Image->YSize;y+=2)        
  234.        for(z=0;z<Image->ZSize;z+=2)
  235.         {
  236.          Sum=0;
  237.          for(i=0;i<2;i++)          
  238.          for(j=0;j<2;j++) 
  239.          for(k=0;k<2;k++) Sum+=Image->Map[x+i][y+j][z+k];
  240.          FindPool_SampleImg3D->Map[x>>1][y>>1][z>>1]=Sum>>3;
  241.         }
  242.       }
  243.      
  244.      Pool->SampledDomainBitmap3D=FindPool_SampleImg3D;
  245.  
  246.      vprintf(stderr,"\n   Domain ");
  247.      Pool->DomainNodes=FindNodes3D(FindPool_SampleImg3D,Pool,B<<1,delta,
  248.                                    FeatureDimensions,(float)TShadeEdge,TRUE);
  249.      
  250.      Pool->DomainSpace=FindSpace3D(Pool,FeatureDimensions,GridRes,(float)TShadeEdge);
  251.  
  252.      return Pool;
  253.     }
  254.    else return (PoolStructure *)NULL;
  255.   }
  256.